      subroutine poisson_fft(ncells, ijkcell, nvtxkm, nvtx, ijkvtx, itdim, iwid, &
         jwid, kwid, precon, tiny, relax, ibp1, jbp1, kbp1, cdlt, sdlt, strait&
         , dx,dy,dz, c1x, c2x, c3x, c4x, c5x, c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y&
         , c6y, c7y, c8y, c1z, c2z, c3z, c4z, c5z, c6z, c7z, c8z, vol, vvol, &
         vvolaxis, q, qtilde, aqtilde, ht, itmax, error, srce, rdt, scalsq, wk&
         , phibc, ex, ey, ez, dive, residu, aq, phi, diag, wkl, wkr, wkf, wke, &
         periodicx)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
      use modify_com_M
      use boundary
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02
!...Switches: -yf -x1
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer  :: ncells
      integer  :: nvtxkm
      integer  :: nvtx
      integer  :: itdim
      integer  :: iwid

      integer  :: jwid
      integer  :: kwid
      integer  :: ibp1
      integer  :: jbp1
      integer  :: kbp1
      integer , intent(in) :: itmax
      real(double)  :: tiny
      real(double) , intent(in) :: relax
      real(double)  :: cdlt
      real(double)  :: sdlt
      real(double)  :: strait
      real(double)  :: dx,dy,dz
      real(double) , intent(in) :: error
      real(double)  :: rdt
      real(double)  :: scalsq
      real(double)  :: wk
      real(double)  :: phibc
      real(double)  :: wkl
      real(double)  :: wkr
      real(double)  :: wkf
      real(double)  :: wke
      logical , intent(in) :: precon
      logical  :: periodicx
      integer  :: ijkcell(*)
      integer  :: ijkvtx(*)
      real(double)  :: c1x(*)
      real(double)  :: c2x(*)
      real(double)  :: c3x(*)
      real(double)  :: c4x(*)
      real(double)  :: c5x(*)
      real(double)  :: c6x(*)
      real(double)  :: c7x(*)
      real(double)  :: c8x(*)
      real(double)  :: c1y(*)
      real(double)  :: c2y(*)
      real(double)  :: c3y(*)
      real(double)  :: c4y(*)
      real(double)  :: c5y(*)
      real(double)  :: c6y(*)
      real(double)  :: c7y(*)
      real(double)  :: c8y(*)
      real(double)  :: c1z(*)
      real(double)  :: c2z(*)
      real(double)  :: c3z(*)
      real(double)  :: c4z(*)
      real(double)  :: c5z(*)
      real(double)  :: c6z(*)
      real(double)  :: c7z(*)
      real(double)  :: c8z(*)
      real(double)  :: vol(*)
      real(double)  :: vvol(*)
      real(double)  :: vvolaxis(*)
      real(double) , intent(inout) :: q(*)
      real(double)  :: qtilde(*)
      real(double)  :: aqtilde(*)
      real(double)  :: ht(*)
      real(double)  :: ex(*)
      real(double)  :: ey(*)
      real(double)  :: ez(*)
      real(double)  :: dive(*)
      real(double)  :: residu(*)
      real(double) , intent(inout) :: aq(*)
      real(double)  :: diag(*)

      real(double), dimension(0:ibp1,0:jbp1,0:kbp1)  :: phi,srce

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	real(double), dimension((ibp1-1)*10) :: wx
	real(double), dimension((jbp1-1)*10) :: wy
	real(double), dimension((kbp1-1)*10) :: wz
	integer :: bcx,bcy,bcz,nx,ny,nz,i,j,l,ijk

nx=ibp1-1
ny=jbp1-1
nz=kbp1-1

! 0 = periodic
! 1 = dirichelet
! other = neumann

if (periodicx) then
bcx=0
else
if(bcMx.eq.1.or.bcMx.eq.0.or.bcMx.eq.-1.or.bcMx.eq.-2) then
 ! dirichlet (sine transform here)
bcx=1
else
 ! neumann
bcx=2
end if
end if

if (periodicy) then
bcx=0
else
if(bcMy.eq.1.or.bcMy.eq.0.or.bcMy.eq.-1.or.bcMy.eq.-2) then
 ! dirichlet (sine transform here)
bcy=1
else
 ! neumann
bcy=2
end if
end if

if (periodicz) then
bcz=0
else
if(bcMz.eq.1.or.bcMz.eq.0) then
 ! dirichlet
bcz=1
else
 ! neumann
bcz=2
end if
end if

write(*,*) 'inizia poisson'

call iniz_poisson(nx,ny,nz,iwid,jwid,kwid,bcx,bcy,bcz,wx,wy,wz)

call poissoneq(nx,ny,nz,bcx,bcy,bcz,dx,dy,dz,srce(1:nx,1:ny,1:nz),phi(1:nx,1:ny,1:nz),wx,wy,wz,rdt, scalsq)

write(*,*) 'fine poisson'

!return

! bc in x
if (bcx==0) then
do j=0,ny+1
do l=0,nz+1
phi(0,j,l)=phi(nx,j,l)
phi(nx+1,j,l)=phi(1,j,l)
enddo
enddo
elseif (bcx==1) then
do j=0,ny+1
do l=0,nz+1
phi(0,j,l)=0.d0-phi(1,j,l)
phi(nx+1,j,l)=0.d0-phi(nx,j,l)
enddo
enddo
else
do j=0,ny+1
do l=0,nz+1
phi(0,j,l)=phi(1,j,l)
phi(nx+1,j,l)=phi(nx,j,l)
enddo
enddo
end if

! bc in y
if (bcy==0) then
do i=0,nx+1
do l=0,nz+1
phi(i,0,l)=phi(i,ny,l)
phi(i,ny+1,l)=phi(i,1,l)
enddo
enddo
elseif (bcy==1) then
do i=0,nx+1
do l=0,nz+1
phi(i,0,l)=0.d0-phi(i,1,l)
phi(i,ny+1,l)=0.d0-phi(i,ny,l)
enddo
enddo
else
do i=0,nx+1
do l=0,nz+1
phi(i,0,l)=phi(i,1,l)
phi(i,ny+1,l)=phi(i,ny,l)
enddo
enddo
end if

! bc in z
if (bcz==0) then
do j=0,ny+1
do i=0,nx+1
phi(i,j,0)=phi(i,j,nz)
phi(i,j,nz+1)=phi(i,j,1)
enddo
enddo
elseif (bcz==1) then
do j=0,ny+1
do i=0,nx+1
phi(i,j,0)=0.d0-phi(i,j,1)
phi(i,j,nz+1)=0.d0-phi(i,j,nz)
enddo
enddo
else
do j=0,ny+1
do i=0,nx+1
phi(i,j,0)=phi(i,j,1)
phi(i,j,nz+1)=phi(i,j,nz)
enddo
enddo
end if

end subroutine poisson_fft


subroutine iniz_poisson(nx,ny,nz,iwid,jwid,kwid,bcx,bcy,bcz,wx,wy,wz)
implicit none
integer :: nx,ny,nz,i,j,l,bcx,bcy,bcz,iwid,jwid,kwid,ijk
real(8), dimension(10*nx) :: wx
real(8), dimension(10*ny) :: wy
real(8), dimension(10*nz) :: wz
real(8) :: dx,dy,dz


! initialization of fft

call trasf_init(nx,wx,bcx)
call trasf_init(ny,wy,bcy)
call trasf_init(nz,wz,bcz)





end subroutine iniz_poisson


subroutine defk(n,bc,k,factor)
implicit none
real(8) :: factor
real(8), dimension(*) :: k
integer :: n,bc

integer ::nn,i
real(8) :: pi

pi=acos(-1.d0)

! bc = 0  => periodic
! bc = 1  => Dirichelet (0)
! bc = 2  => Neumann (0)
! in all cases solves only for unknown.
! periodic -> puts periodicit to the first outside
! dirichelet -> puts to zero the point outside to the left and the the right
! neunmann -> sets derivative to zero on first and last point in the array

if(bc==0) then
! fft standard
if(mod(n,2)==0) then
nn=n/2
else
nn=(n+1)/2
end if

k(1)=0.d0
do i=2,nn
k(2*i-2)=(i-1)*2.d0*pi/dble(n)
k(2*i-1)=(i-1)*2.d0*pi/dble(n)
enddo

if(mod(n,2)==0) then
k(n)=1.d0
end if

factor=1.d0/dble(n)


else if (bc == 1) then
! sine fft

do i=1,n
k(i)=i*pi/dble(n+1)
enddo

factor=1/2.d0/(dble(n+1))

else
! cosine fft

do i=1,n
k(i)=(i-1)*pi/dble(n-1)
enddo

factor=1/2.d0/(dble(n-1))

end if

end subroutine defk

subroutine trasf_init(n,w,bc)
implicit none
real(8), dimension(*) :: w
integer :: n,bc

if(bc == 0) then
CALL DFFTI (N,W)
else if(bc==1) then
call dsinti(n,w)
else if(bc==2) then
call dcosti(n,w)
end if

end subroutine trasf_init

subroutine trasf(n,bc,vec,w)
implicit none
integer :: n,bc
real(8), dimension(n) :: vec
real(8), dimension(*) :: w

if (bc == 0) then
CALL DFFTF (n,vec,w)
else if (bc == 1) then
call dsint (n,vec,w)
else
call dcost (n,vec,w)
end if

end subroutine trasf


subroutine trasf_inv(n,bc,vec,w)
implicit none
integer :: n,bc
real(8), dimension(n) :: vec
real(8), dimension(*) :: w

if (bc==0) then
call dfftb(n ,vec,w)
else if (bc == 1) then
call dsint (n,vec,w)
else
call dcost (n,vec,w)
end if

end subroutine trasf_inv


subroutine poissoneq(nx,ny,nz,bcx,bcy,bcz,dx,dy,dz,rho,phi,wx,wy,wz,rdt, scalsq)
implicit none
integer :: nx,ny,nz,bcx,bcy,bcz
real(8), dimension(nx,ny,nz) :: rho,phi
real(8), dimension(10*nx) :: wx
real(8), dimension(10*ny) :: wy
real(8), dimension(10*nz) :: wz
real(8) :: dx,dy,dz,rdt, scalsq

integer :: i,j,l
real(8) :: factorx,factory,factorz,delta2
real(8), dimension(nx) ::  kx
real(8), dimension(ny) ::  ky
real(8), dimension(nz) ::  kz


call defk(nx,bcx,kx,factorx)
call defk(ny,bcy,ky,factory)
call defk(nz,bcz,kz,factorz)

phi=rho

do j=1,ny
do l=1,nz
CALL trasf(nx,bcx,phi(:,j,l),wx)
enddo
do i=1,nx
CALL trasf(nz,bcz,phi(i,j,:),wz)
enddo
enddo
do i=1,nx
do l=1,nz
CALL trasf(ny,bcy,phi(i,:,l),wy)
enddo
enddo

do i=1,nx
do j=1,ny
do l=1,nz
delta2=-2.d0*((cos(kx(i))-1.d0+1.d-5)/dx**2+ &
              (cos(ky(j))-1.d0+1.d-5)/dy**2+ &
              (cos(kz(l))-1.d0+1.d-5)/dz**2)
delta2=-delta2*scalsq-rdt
!delta2=4.d0/dx**2*(cos(ky(j)/2.d0)*cos(kz(l)/2.d0)*sin(kx(i)/2.d0))**2 + &
!       4.d0/dy**2*(cos(kx(i)/2.d0)*cos(kz(l)/2.d0)*sin(ky(j)/2.d0))**2 + &
!       4.d0/dz**2*(cos(kx(i)/2.d0)*cos(ky(j)/2.d0)*sin(kz(l)/2.d0))**2

phi(i,j,l)=phi(i,j,l)/delta2
enddo
enddo
enddo

do j=1,ny
do l=1,nz
CALL trasf_inv(nx,bcx,phi(:,j,l),wx)
enddo
do i=1,nx
CALL trasf_inv(nz,bcz,phi(i,j,:),wz)
enddo
enddo
do i=1,nx
do l=1,nz
CALL trasf_inv(ny,bcy,phi(i,:,l),wy)
enddo
enddo

phi=phi*factorx*factory*factorz

end subroutine poissoneq
